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Abstract 

We discuss lattice simulations of light nuclei at leading order in chiral effective field 
theory. Using lattice pion fields and auxiliary fields, we include the physics of instantaneous 
one-pion exchange and the leading-order S-wave contact interactions. We also consider 
higher-derivative contact interactions which adjust the S-wave scattering amplitude at higher 
momenta. By construction our lattice path integral is positive definite in the limit of exact 
Wigner £77(4) symmetry for any even number of nucleons. This ££/(4) positivity and 
the approximate SU(A) symmetry of the low-energy interactions play an important role 
in suppressing sign and phase oscillations in Monte Carlo simulations. We assess the 
computational scaling of the lattice algorithm for light nuclei with up to eight nucleons and 
analyze in detail calculations of the deuteron, triton, and helium-4. 



I. INTRODUCTION 

The underlying theory of strong interactions, quantum chromodynamics (QCD), describes 
the interactions of quarks and gluons. While analytic calculations of the properties of 
confined quarks and gluons inside hadrons are not possible, a model-independent way of 
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calculating observables directly from QCD is provided by lattice field theory. Recent ad- 
vances in lattice QCD have made it possible to calculate the spectrum and properties of 
various isolated hadrons. There has also been progress in calculating low-energy hadronic 
interactions such as pion-pion scattering Other hadronic interactions such as 

nucleon-nucleon scattering are more difficult, but there has been some promising recent work 
in this direction as well p, [(], Q] • 

Unfortunately lattice QCD calculations of many-body systems of nuclear and neutron 
matter or even few-body systems beyond two nucleons are presently out of reach. Such 
simulations would require pion masses at or near the physical mass and lattices several times 
longer in each dimension than used in current simulations. But the greatest challenge would 
be to overcome the exponentially small signal-to-noise ratio for simulations at large quark 
number. For many-body systems this is manifested as complex phase oscillations when 
adding a quark chemical potential. For few-body systems the calculation can be done at zero 
chemical potential by measuring correlation functions involving 3A-quark operators, where 
A is the number of nucleons. However here the signal-to-noise problem reappears in the 
antisymmetrization over quarks and in the small overlap between Monte Carlo configurations 
for the QCD vacuum versus the A-nucleon ground state. 

For few- and many-body systems in low-energy nuclear physics one can make further 
progress by working directly with hadronic degrees of freedom. There are several possible 
choices for the form of the nuclear forces and the computational methods used to describe 
the interactions of low-energy protons and neutrons. 

For systems with four or fewer nucleons, a numerically exact approach is provided by 
the Faddeev-Yakubovsky integral equations. Three-nucleon continuum observables as well 



as the triton and a-particle binding energies 



were extensively studied within the Faddeev- 



Yakubovsky scheme based on a variety of modern semi-phenomenological nucleon-nucleon 
potential models including the CD-Bonn [9], CD-Bonn 2000 [h]], Argonne V18 U| and 
Nijmegen |l2| potentials. Three-nucleon forces were also incorporated using the Tucson- 
Melbourne [3, [jjj], Urbana-IX [ijj and other models. For a comprehensive review on the 



calculations in the three-nucleon continuum the reader is referred to [16]. The same compu- 
tational scheme was applied to nuclear forces derived in chiral effective field theory (ChEFT) 
both at next-to- leading order (NLO) 17| and at next-to- next-to-leading order (NNLO) 18] 
in the chiral expansion. Applications of the low-momentum interaction potential V[ ow k 
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2 11 ] to few- nucleoli systems are considered in Refs. [22j, |23| . Further computational 



techniques such as, e.g., the expansion in hyperspherical harmonics 



24|, the Lorentz inte- 



gral transform method 251, the stochastic variational method 26] and the Kohn- variational 

n 

approach [27( were also applied to few-nucleon systems. 

For systems with more nucleons one must rely on techniques such as Monte Carlo simula- 
tions or basis-truncated eigenvector methods. There have been a number of Green's Func- 



tion Monte Carlo simulations o 
potentials, see for example [15 



ight nuclei based on AV18 as well as other phenomenological 
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331 ] . A related technique implementing 



diffusion Monte Carlo with auxiliary 



tron matter and neutron droplets 



Mds has been used to study the ground state of neu- 
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331 • The No-Core Shell Model (NCSM) is a 



different approach to light nuclei which uses basis-truncated eigenvector methods. There 



have been several N 
models, cf. 



38, 
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CSM calculations using various different phenomenological potential 
rere are also NCSM calculations which have used nuclear 
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forces derived from ChEFT [42j, |43|] . Quite recently there has been work in constructing a 
low-energy effective field theory within the framework of truncated basis states used in the 
NCSM formalism 44|. A benchmark comparison of many of the methods listed above as 



451 ] . A review article on various methods used for 



well as other techniques can be found in 
light nuclei can be found in 461 ] . 

In this study we consider nuclear lattice simulations of light nuclei using chiral effective 
theory. The nuclear lattice approach addresses the few- and many-body problem in nuclear 
physics by applying non-perturbative lattice methods to low-energy nucleons and pions. 
The chiral effective Lagrangian is formulated on a Euclidean lattice and the path integral is 
evaluated by Monte Carlo sampling. Pions and nucleons are treated as point-like particles 
on the lattice sites, and tc times the inverse lattice spacing sets the cutoff scale in momentum 
space. By using hadronic degrees of freedom and concentrating on low-energy physics, it is 
possible to probe larger volumes, lower temperatures, and far greater numbers of nucleons 
than in lattice QCD. In some cases the sign and complex phase oscillations in Monte Carlo 
simulations can be significantly reduced or even completely eliminated. 

The first study combining lattice methods with effective field theory for low-energy nuclear 
ohysics looked at infinite nuclear and neutron matter at nonzero density and temperature 



471 ] . The approach we use here is based on chiral effective field theory starting at leading 



order. This lattice formalism was also used in [48] to study neutron matter at nonzero 
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temperature. We list some features of the nuclear lattice approach which seem promising 
and distinguish it from other few- and many-body techniques. 

One unique feature of the lattice effective field theory approach is the ability to study 
in the same formalism both few- and many-body systems as well as zero- and nonzero- 
temperature phenomena. A large portion of the nuclear phase diagram can be studied 
using exactly the same lattice action with exactly the same operator coefficients. A second 
feature is the computational advantage of many efficient Euclidean lattice methods developed 
for lattice QCD and condensed matter applications. This includes the use of Markov 
Chain Monte Carlo techniques, Hubbard- Stratonovich transformations 
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501 ] . and non- 



local updating schemes such as a hybrid Monte Carlo (5JJ. A third feature is the close 
theoretical link between nuclear lattice simulations and chiral effective field theory. One can 
write down the lattice Feynman rules and calculate lattice Feynman diagrams using precisely 
the same action used in the non-perturbative simulation. Since the lattice formalism is 
based on chiral effective field theory we have a systematic power-counting expansion, an a 
priori estimate of errors for low-energy scattering, and a clear theoretical connection to the 
underlying symmetries of QCD. 

Nuclear lattice simulations were used to study the triton at leading-order in pionless 
effective field theory with three-nucleon interactions [52[]. In the present investigation we 
consider the physics of instantaneous one-pion exchange and the leading-order S-wave con- 
tact interactions. We also consider higher- derivative contact interactions which adjust the 
S-wave scattering amplitude at higher momenta. We calculate binding energies, radii, and 
density correlations for the deuteron, triton, and helium-4, and probe the computational 
scaling in systems with up to eight nucleons. 

II. CHIRAL EFFECTIVE FIELD THEORY IN THE FEW-NUCLEON SECTOR 

Chiral perturbation theory in the purely mesonic sector has a rigorous chiral counting 



scheme. In the one-nucleon sector a chiral counting scheme can be established 



means such as the heavy-baryon formulation 
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54[ | or infrared regularization 55j . In each 



Dy various 



case Green's functions are expanded in increasing powers of pion masses and small momenta, 
and the chiral expansion corresponds to a loop expansion. 

In the few-nucleon sector however one has to deal with non-perturbative effects. Pertur- 
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bation theory fails at low energies in the few-nucleon sector due to enhanced contributions 
from reducible diagrams which contain purely nucleonic intermediate states. In order to 
circumvent this problem, one derives first the interaction kernel (or effective potential) from 



all possible irreducible terms without purely nucleonic intermediate states [56j, [57|. The 
interaction kernel does not contain small energy denominators and obeys the conventional 
chiral counting scheme. The Green's function is then obtained by iterating the interaction 
kernel to infinite order in a bound state or scattering state equation. 

At lowest order in the chiral expansion the effective Lagrange density is 

11 V 2 
C = -d u 7T ■ d^ir - -mlir 2 + NHd N + iV f — JV 
2 2 n ° 2m 

- -^-iVV a ■ -VirN - -C{N^N)(N^N) - -Cj{N^tN) ■ (iWiV). (1) 
2/7T 2 2 



We use the same notation as used in [581 ] . N is the nucleon field with spin and isospin 
degrees of freedom. The vector arrow in a signifies the three-vector index for spin. The 
boldface for r and w signifies the three-vector index for isospin. We take for our physical 
constants m = 938.92 MeV as the nucleon mass, m n = 138.08 MeV as the pion mass, 
/„- = 93 MeV as the pion decay constant, and gA = 1-26 as the nucleon axial charge. In 
order to reduce sign and complex phase oscillations in the Monte Carlo calculation with 
auxiliary fields (see 59J]) we work with the leading-order contact interactions C and Cj 



rather than the more standard interaction coefficients Cs and CV corresponding with 

- ^C S (N^N)(N^N) - —Ct(N^(tN) ■ (N^aN). (2) 
Both forms for the interactions are exactly the same if we set 

C = Cs — 2C7 1 , Cj = — Ct- (3) 
From the effective Lagrangian in ([T]) the NN effective potential can be derived by appl ying 



the method of unitary transformations 60j, Q-box expansion 6l|, or other techniques 62, 



631 ] . At leading order (LO) the NN effective potential consists of the two contact interactions 



and instantaneous one-pion exchange, 



V LO = C + C I r 1 .r 2 -[—j r x -r 2 - - - , (4) 
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where q = p ' — p is the nucleon momentum transfer. We can reproduce the desired iteration 
of Vlo if we start with the Lagrange density, 

11 V 2 

C = — Vtt • -Vtt - -mln 2 + NUd N + N j N 

2 2 2m 

- ^-iVV a ■ -VizN - -C(N^N)(N^N) - -C^N^tN) ■ (NWN), (5) 
2/tt 2 2 

and evaluate the NN scattering amplitude non-perturbatively. We note that the pions have 
no time derivatives. Therefore they can only be exchanged instantaneously between nucleons 
and do not propagate in time. Clearly the Lagrangian in Eq. (J5D is not valid for external pion 
fields. The two-nucleon Green's function derived from the path integral representation with 
this Lagrangian reproduces the solution of the corresponding Lippmann-Schwinger equation 
with the leading order effective potential. Another advantage of treating pions this way is 
that the nucleon self-energy exactly vanishes and the nucleon mass is not renormalized. 



III. LATTICE NOTATION 



In this study we assume exact isospin symmetry and neglect electromagnetic interactions. 
We use n to represent integer- valued coordinate lattice vectors on a 3 + 1 dimensional space- 
time lattice and k to represent integer- valued momentum lattice vectors. A subscripted V 
such as in n s represents purely spatial lattice vectors. denotes the unit lattice vector in 
the time direction, and l s — 1, 2, 3 are unit lattice vectors in the spatial directions, a is 
the spatial lattice spacing, L is the length of the cubic spatial lattice in each direction, a t is 
the lattice spacing in the temporal direction, and L t is the length in the temporal direction. 
We define a t as the ratio between lattice spacings, a t = a t j 'a, and define h = a t /(2m). 
Throughout we use dimensionless parameters and operators, which correspond with physical 
values multiplied by the appropriate power of a. Final results are presented in physical units 
with the corresponding unit stated explicitly. 

To avoid confusion we make explicit in our lattice notation all spin and isospin indices. 
We use c and c* to denote the anticommuting Grassmann variables for the nucleons and a 
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and a' to denote annihilation and creation operators. We use the subscript notation 







c 0,0 




a 1,p 




a 0,0 


c i,p 




Cl,0 




a l,p 




a l,0 






Co.l 








a>o,i 


. C Un . 




. c l,l . 




. a Un . 




. a i,i . 



(6) 



The first subscript is for spin and the second subscript is for isospin. We use 17 with 
/ = 1,2,3 to represent Pauli matrices acting in isospin space and as with 5 = 1,2,3 to 
represent Pauli matrices acting in spin space. We note that on the lattice the spin symmetry 
is reduced to the cubic subgroup 50(3, Z) of 50(3) while isospin symmetry remains intact 
as the full SU (2) symmetry. 

IV. PATH INTEGRAL FOR FREE NUCLEONS AND INSTANTANEOUS PIONS 

If we could take the continuum limit, the accuracy of our calculation would be determined 
by the order k we choose to truncate the chiral expansion. This would correspond with 
(p/A x ) h , where p is a typical low-energy momentum scale and A x = 47tf n ~ 1.2 GeV the 
scale of spontaneous chiral symmetry breakdown. However taking the continuum limit is 
not possible due to the non-perturbative treatment of the chiral effective Lagrangian on the 
lattice as this would require an infinite set of counterterms. The lattice cutoff A must be 
chosen to remain below the scale A x . This in turn introduces an error of the general form 
(p/A) fcl (A/A x ) fe2 due to the finite cutoff and missing counterterms. Even for the free nucleon 
case finite cutoff errors occur which can be traced back to the discretized lattice propagator. 
In this analysis we use an 0(a 4 )-improved action for the lattice kinetic energy, as shown in 
Figure [H 

Similarly, the interactions in the continuum limit can be organized in the chiral expansion 
as leading order, next-to-leading order, etc. But again for a chosen chiral order we may wish 
to include additional improvements to the interactions which reduce the finite cutoff errors. 
In this analysis we start by considering the simple LO lattice action without improvement, as 
shown in Figure [2j We note that the diagram shown is a bit simplistic since the improvement 
terms may in general include corrections for effects at nonzero lattice spacing such as broken 
Galilean invariance, etc. 
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FIG. 1: Suppression of finite cutoff errors by introducing improved actions for the nucleon kinetic 
energy. 



Throughout our discussion we consider both the path integral formalism and the transfer 
matrix formalism. The path integral formalism is useful for deriving the lattice Feynman 
rules and auxiliary field formulations, while the transfer matrix is used for the Monte Carlo 
simulations of light nuclei. We start with the path integral formalism. Let Z$ N be the 
lattice partition function for free nucleons 




(7) 



where 



DcDc* 



Y[ dc i,j(n)dclj(n). 



(8) 
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FIG. 2: Hierarchy of interactions according to chiral order as well as improvements to each order 
due to nonzero lattice spacing. 

We use an 0(a 4 )-improved lattice action with next-to- next-to-nearest neighbor hopping, 



We expect the 0(a 4 )-improvement in the lattice dispersion relation to be useful when mea- 
suring scattering phase shifts on the lattice. 

In this leading-order study we consider instantaneous one-pion exchange and no other 
interactions involving pions. In our lattice formalism the pion field does not propagate 
in time and does not couple to physical pions. This allows to avoid the problem of non- 
perturbative dynamical pion fields producing unrenormalized pion loops to all orders. If 
at some point later on we wish to include interactions with physical low-energy pions we 




ri,l s ,i,j 



oh r 

+ [c^(ff)c iii (n + 2Z a )+c^(fr)c iii (fr-2Z.) 



n,l a ,i,j 



(9) 



simply insert the corresponding operators with external pion fields. 

The lattice action for free pions with purely instantaneous propagation is 

Smim) = at(r%- + 3) E ^2,K I (n)'K I (n) - a t E ^^A^Ti^n + l s ). (10) 

7=1,2,3 n 7=1,2,3 n,l s 

In order to simplify the Monte Carlo updating scheme later in our discussion it is helpful at 
this point to define a rescaled pion field, 



(11) 



where 



Then 



QiT = a t (ml + 6). 



7=1,2,3 

In momentum space the action is 



s,M) = -^- 3 E £^(-*K(*) 



7=1,2,3 n,l s 



1 a t 



L t L 3 

7=1,2,3 % 

and so 

|L>7r> 7 (^X(0)exp [S^] 



2wk, 



UL t V - 

^ / cos ■ T 



f Dvr^exp [-S wn ] 
where the pion propagator is 



no sum on 7) 



E 



1-— E cos i , 

i s =l,2,3 



27rfc, 



The pion correlation function at spatial separation n s is then 



, ... , / £>7r>£ (ra,X(5) exp [S n7T ] 
7rj(n s )7rj(0) > = t^— - (no sum on 7) 



J TJtTj exp [— £?„ 



L 3 ^ 



e l - 



(12) 



(13) 



(14) 



e-^^e-^^DM, (15) 



(16) 



(17) 
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V. PION-NUCLEON COUPLING 



There are various ways to introduce spatial derivatives of the pion field on the lattice. The 
simplest definition for the gradient of tt'j is to define a forward-backward lattice derivative. 
For example we can write 

d 1 n' I (n) = ±[n' I (n + i)-n' I (n-l)]. (18) 



This is the method used in 



48 



]. The disadvantage is that it is a coarse derivative involving 
a separation distance of two lattice units. We can avoid this if we think of the pion lattice 
points as being shifted by —1/2 lattice unit from the nucleon lattice points in each of the 
three spatial directions. For each nucleon lattice point n n ucieon we associate a pion lattice 
point n pion , 



1* 1 4 1 



n 



pion 



^nucleon „1 „2 3. (19) 



2 2 2 



Then we have eight pion lattice points forming a cube centered at n nuc i e0 n, 

^pion; ^pion + 1; ^-pion + 2, ^pion + 3, 

^pion + 1 + 2, n pion + 2 + 3, n pion + 3 + 1, n pion + 1 + 2 + 3. (20) 

We use the same lattice vector notation n for both nucleons and pions. However for nucleon 
fields and auxiliary fields to be introduced later n represents n n ucicon while for pion fields n 
refers to n P i on . 

The eight vertices of the pion cube in (120]) can be used to define spatial derivatives of the 
pion field. For each spatial direction S — 1, 2, 3 we have 

A s n' I (n) = ^ J2 (-~L) Us+ %(n + v), * = v{L + u 2 2 + u 3 3. (21) 

^1,^2,^3=0,1 

The lattice pion-nucleon coupling in our lattice action is 

S^ N {^c,c*) = ^= ^sAWpsM, (22) 

S,/=l,2,3 

where ps,i(n) is the spin-isospin density, 

Ps,i(n)= Yl c ij(^) \. a s\u' [ T i]jf Ci'j'in)- (23) 
«,i,j',i'=o,i 
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VI. S-WAVE CONTACT INTERACTIONS 



There are two S-wave contact interactions at lowest order. Following 59] we choose the 
form 



S 



NNNNK^ 



C, C 



(24) 



ra 7=1,2,3 n 

where p(n) and p/(n) are the S77(4)-symmetric and isospin densities respectively, 



P(n) = Y c lj^) c 

»J'=0,1 



(25) 
(26) 



Since the isospin singlet channel is more strongly attractive than the isospin triplet channel, 
we anticipate the signs for these coefficients to be C < and Gj > 0. This will be confirmed 
in Sec. IVIIIt where the two-nucleon system is studied in detail. As noted above, these can 
be written in terms of the more familiar coefficients Cs and Ct using the identity 



C — C,s — 2Ct, 



Cj — — Cj'. 



(27) 



We can use the Gaussian integral identities 

fOO 



<>X P 1 —7T [Pi™)? 



1 



'2tt 



dsexp 



--s 2 + y/-Ca t p{n) ■ s 



and 



exp 



C T a t 



Y WW*) 

1=1,2,3 ) 



I eX P 



\7=1,2,3 



271 



-\ Y (si? + iVCi<*t Y 



Pi(n) ■ S! 



7=1,2,3 



7=1,2,3 



Let us define the auxiliary field actions, 

s ss (s,s/) = ^Y s2 ^ + - y Y^ Sl ^ 



2 t-^i K ' 2 

n 7=1,2,3 n 



S sNn{ s , g 7, c, c*) = -y/-Ca t ^ s(n)p(n) - i\j Cia t Y /^ij^Pii^)- 

n 7=1,2,3 n 

Then we have 

y DsDs T exp [-^(s, s/) - S s n N (s, s 7 , c, c*)] oc exp [-Sf} N fi N (c, c*)] , 



(28) 



(29) 



(30) 
(31) 

(32) 
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where 



n . 



DsDsj = Y\_ds(n)dsi( 

n,I 

If we put all the pieces together the full path integral action at leading order is 
Zlo oc J DcDc* Dir'fDsDsi exp [— Slo(c, c*, tx\, s, s/)] , 



(33) 



where 



Slo — $nn + S n7T + S wNN + S ss + S sNN . 



(34) 



(35) 



VII. TRANSFER MATRIX WITH AUXILIARY FIELDS 

The transfer matrix is the analog at nonzero temporal lattice spacing of the operator 
exp (—HAt). In order to derive the transfer matrix corresponding with the path integral 
action Slo we use the correspondence [64j, [65] 



Tr | : F Lt -i a\, r {n' s ), a itj {n s )\^ : x • • • x : F al tj ,(n' s ), a it j{n s ) 
= / DcDc*exp c iA™*> n t) [ c i,j(^s,n t ) - Cij{n s ,n t + 1)] \ 

I nt=0n a ,i,j \ 
L t -1 

x n 



(36) 



n t =0 



for general functions Fi and antiperiodic boundary conditions in the time direction, 
Cij(n s ,L t ) = — Ci j(n s , 0). The : : symbols in ( |36|) denote normal ordering. Let us de- 
fine the S'L r (4)-symmetric, isospin, and spin-isospin densities written in terms of creation 
and annihilation operators, 



P a ' a (n. 



M=0,1 



Pi K n s) ^ 

i,j,f =0,1 



T i]jj> a i,j'(n s ), 



(37) 

(38) 
(39) 



i,j,i',j'=0,l 
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We also define the 0( a 4 ) -improved free nucleon lattice Hamiltonian 

49 



H ivcc = — E a \,j{^s)a%,j{n s ) 



Am 



E a lj(^s)cii,j{n s + l s ) + al j (n s )a id (n s - l s ) 



+ E [ a L-("s) a *j'(^> + + a L("s)aij(n s - 2l s ) 



,1 



— £ 

180m ^ 



a\ j {n s )a id {n s + 31 s ) + a\j(n s )a itj {n s - 3l s 



(40) 



Using (|36j) the path integral with auxiliary fields can be expressed in the transfer matrix 
formalism as 



Z LO oc j DitjDsDsj exp [S W7T - S ss ] x Tr {M {Lt ~ l) (n'j, s, sj) x • • • x M (0) «, s, Sj )} 



(41) 



where 



M^\tt' i ,s,s i )= : expJ-^ free a t -^|=^A S 7rKn s ,n t )P5!/ a (^) 

I s,i 

+ ^-Ca t J2 s@» nt)p a '' a (n s ) + ^^/C^~ t J2Yl n t)pf' a (n s ) } ( 12) 



VIII. THE TWO-NUCLEON SYSTEM 



For the two-nucleon system the entire linear space is small enough for typical lattice 



re lattice using iterative sparse 



66j |. To do this calculation we 



volumes that we can find the low-energy eigenstates on t 
matrix eigenvector methods such as the Lanczos method 
construct the transfer matrix with only nucleon fields. It is convenient to define 

/ D^AsM^AsMO) exp [S v 



G Sl s 2 (n a ) 



(no sum on /) 



/ Dix\ exp 

= y Q E E (-irM-i)^(^ + ^fW(6) 

The path integral can now be written as 



(43) 



Z LO oc Tr {M {Lt - 1] x • • • x 
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M®} 



(44) 



where 



= : exp \ -H ircc a t - l -Ca t ^ \p a ^ a (n s )\ * - ±da t £ £ [pf ' a 



J ft. 



+ 1^ G s 1 s 2 (n s ,i -n s ,2)pi;}{n a> x)pi;}{n Sj2 ) } : . (45) 

Si,S2,I n s> i,n Si2 

We now calculate the two-nucleon spectrum in a periodic cube of length L and use this 
information to determine the contact interaction coefficients C and Cj. We make use of 
0, H, 68 1 which relates the two-particle energy levels in a periodic cube 



Liischer's formula 



of length L to the S-wave phase shift, 



pcot5 (p) = —S (rj) 



Lp 

"=1 to 



where S(rj) is the three-dimensional zeta function, 

6(A 2 - n 2 ) 



S(rj) = lim 

A— >oo 



E 



47rA 



rr — rj 



(46) 



(47) 



For |r?| < 1 we can expand in powers of rj 

1 



lim 

Tj A^oo 



E 



fl(A 2 - n 2 
n 2 — 77 



- 47rA 



- + S + -W + S 2 77 2 + S377 3 • • • 
V 



where 



Sr, = lim 



A^oo 



E 



#(A 2 - n 2 ^ 



rr 



47rA 



5j ^ (n 2 V +1 



i>i- 



(48) 
(49) 

(50) 
(51) 



The first few coefficients are 



S = -8.913631, St = 16.532288, S 2 = 8.401924, S 3 = 6.945808, 

S 4 = 6.426119, S 5 = 6.202149, S 6 = 6.098184, S 7 = 6.048263. (52) 

Liischer's formula does not include cutoff effects or the contribution from coupled higher 
partial waves for particles with spin. However we can neglect such corrections at asymp- 
totically small momenta. For small momenta we have the effective range expansion, 

1 1 



p cot 5 (p) 



+ 7; r oP + 

Oscatt ^ 



(53) 
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(54) 



where a scat t is the scattering length and r is the effective range. In terms of rj, the energy 
of the two-body scattering state is 

E = f_ = V_ /2tt X ~ 
m m \ L 

S is an analytic function of i] near r] — 0, and so we can consider both E < and E > 0. 
We decouple the spin-singlet and spin-triplet contact interactions by expressing C and Cj 
as a linear combination of coefficients Cis and Cs Sl , 



C = (3Ci So + Cs Sl )/4, 
C I = (Ci So -Cs Sl )/4. 



(55) 
(56) 



The value of Ci Sl is tuned to give the physical deuteron binding energy, —2.224575(9) MeV. 
The value of Ci So is tuned using Luscher's formula to give the physical 1 So scattering length, 
-23.76(1) fm. 

At leading order in the two-nucleon system we expect finite cutoff errors to scale roughly 
as 0(A _1 ) or 0(a). On the lattice we can relate this cutoff error to the probability that 
both nucleons occupy the same lattice site. This localized two-nucleon state has a large 
positive expectation value for the kinetic energy and a large negative expectation value for 
the potential energy. Let _gi° callzcd be the expectation value of the total energy. £ , i° callzed 
need not be small compared with the cutoff energy A 2 / (2m). Therefore transfer matrix 
elements involving this state may have a significant dependence on the temporal lattice 
spacing even for a^ 1 as large as the cutoff energy. This dependence shows up clearly in the 
0(A _1 ) cutoff error, and we see the effect in the following results. 

For a = (100 MeV)- 1 and a t = (70 MeV)" 1 , (200 MeV)- 1 , (10000 MeV)" 1 we set the 
coefficients C , 35 1 and Cis using Luscher's formula. We also use Luscher's formula to deter- 
mine the effective range for the 1 S partial wave and both the scattering length and effective 
range for the 3 S'i partial wave. The results are shown in Table 1. 

Table 1: Coefficients and S-wave parameters for a = (100 MeV) -1 





a Sl (MeV- 2 ) 


a So (MeV- 2 ) 


rl So (fm) 


a s ?att (fm) 


3 Q 

r 1 (fm) 


a t = (70 MeV) -1 


-5.714 x 10- 5 


-5.021 x 10- 5 


-0.179(7) 


4.153(5) 


-0.48(2) 


at = (200 MeV)" 1 


-6.706 x 10- 5 


-5.794 x 10- 5 


0.71(2) 


4.522(1) 


0.30(2) 


at = (10000 MeV)- 1 


-7.151 x 10- 5 


-6.126 x 10- 5 


1.03(2) 


4.664(1) 


0.53(2) 


experiment 






2.75(5) 


5.424(4) 


1.759(5) 
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Note that the values for Cis and C3$ 1 are reasonably close to the ones found at NLO 



and NNLO in the continuum formulation 



691 ] . Also, within the pionless framework both 



values should be identical in the Wigner symmetry limit. The errorbars on the lattice data 
in Table 1 are error estimates from the least squares fit. Since we work at leading order we 
do not expect agreement with the experimental values for the effective ranges. However it 
is interesting to note that the effective ranges are actually negative for the largest temporal 
lattice spacing. In the following we explain how this happens. 

For some small fixed Euclidean time interval At consider all transition amplitudes between 
two-nucleon states. If a t <C At then there are many temporal lattice steps in the time 
interval At. Any transition amplitude involving states with two nucleons close together is 
enhanced to some degree by the negative potential energy of the delta function potential. 
On the other hand if a t = At then there is only one temporal lattice step. In this case 
only the forward matrix element for incoming and outgoing localized two-nucleon states 
is enhanced by the delta function potential. This produces a sharp central peak in the 
two-nucleon wavefunction where the nucleons overlap and explains the decrease in effective 
range. By the same reasoning we also expect a smaller value for the deuteron root-mean- 
square radius r^. In Table 2 we show results for and the deuteron quadrupole moment, 
Qd, along with corresponding experimental values. The experimental value quoted for rj is 
for the point proton root-mean-square radius. 

Table 2: Properties of the deuteron for a = (100 MeV) -1 





r d (fm) 


Q d (fm 2 ) 


at = (70 MeV)- 1 


1.566(1) 


0.144(1) 


at = (200 MeV)" 1 


1.668(1) 


0.171(1) 


at = (10000 MeV)" 1 


1.736(1) 


0.179(1) 


experiment 


1.9671(6) 


0.2859(3) 



As expected the root-mean-square radius of the deuteron is smaller than the physical value, 
and the deviation is greater for larger values of at- The smaller radius also results in a 
substantial reduction in the quadrupole moment. 
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IX. ZERO-RANGE CLUSTERING INSTABILITY 



We have found that the deuteron wavefunction at leading order shows some deficiencies 
which presumably get fixed at higher order in chiral effective field theory. But since we now 
have in hand the coefficients of the leading-order contact interactions for lattice spacing a = 
(100 MeV)- 1 and a t = (70 MeV)~\ (200 MeV)- 1 , (10000 MeV)^ 1 , we can consider systems 
with more than two nucleons at leading order in chiral effective field theory. Unfortunately 
here we find more problems. In the helium-4 system we discover that the ground state is 
severely overbound and consists almost entirely of the quantum state with all four nucleons 
occupying the same lattice site. This clustering instability can be understood as the result 
of two contributing factors. The first is that chiral effective field theory at leading order 
gives a poor description of S-wave scattering above a center of mass momentum of 50 MeV. 
The leading-order contact interactions are momentum independent and, as a result, are 
too strong at high momenta. The second is a combinatorial enhancement of the contact 
interactions when more than two nucleons occupy the same lattice site. This effect has 
been studied in two-dimensional large- N droplets with zero-range attraction [70J. Similar 



effects have a 
lattices 



71 



so been considered in systems of higher-spin fermions in optical traps and 
721 ]. To illustrate we briefly discuss how the problem arises in 5'L r (4)-symmetric 
pionless theory at leading order using a Hamiltonian lattice formalism. 

Let £'J ocallzed be the expectation value for the kinetic energy of a single nucleon localized 
on a single lattice site and let Vi < be the potential energy between two nucleons on the 
same lattice site. If we fix the two-particle scattering length then both £ , j ocallzed anc l y 2 
scale linearly with the cutoff energy, 

A 2 

^localized _ _y A = 7TO _ . (57) 

1 2m' y J 

A detailed calculation of V2 for infinite scattering length can be found in [73[]. The total 
energies associated with putting two, three, or four nucleons on the same lattice site are 

^localized = ^localized + (gg) 
^localized = 3 £localized + 3^ ( 5g ) 
^localized = ^localized + (qq) 

An instability forms as we increase A because the kinetic energy scales as the number of 
nucleons, A, while the potential energy scales as (^) . Of course the Pauli exclusion principle 
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prevents more than four nucleons from sitting on the same lattice site, and so the problem 
is most severe in the four-nucleon system. 



In the leading-order pionless theory it has been shown that V2 < — E l ° callzed |73|. There- 
fore £ji° callzed is negative and scales with the cutoff energy. This produces an instability for 
the three-nucleon system in the absence of three-body forces or other stabilizing effects. The 
instability of the triton for zero range forces was first studied by Thomas [7^]. There have 
been a number of more recent studies of the triton in pionless effective field theory as well as 



more general three- 



52 



75 



76 



77 



78 



jody systems with short range interactions and long scattering lengths 



79 
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8 If . It has also been shown that when the cutoff dependence 



in the three-nucleon system is removed using a three-nucleon contact interaction, then the 



binding energy of the four-nucleon system appears also to be cutoff independent |82|, [83]. 
In our lattice Hamiltonian notation we denote V3 as the potential energy associated with 
the three-nucleon contact interaction. The new localized energies are then 

^localized = ^localized + y^ (g^ 

^localized = g^localizod + gyr + y^ (62) 

^localized = ^localized + gyr + ^ (gg) 

Clearly ^ ocallzed would be stabilized by 4V3 for sufficiently large V3 > 0. However for realistic 
nuclear binding energies, it was found that the desired cutoff independence in helium-4 does 
not emerge until the cutoff momentum A exceeds 8 fm -1 jsjij]. Unfortunately this high cutoff 
momentum makes it a difficult starting point for lattice simulations of realistic light nuclei. 
A cutoff momentum of 8 fm -1 corresponds with a lattice spacing of about 0.4 fm. From 
a computational standpoint this combination of short lattice spacing and strong repulsive 
forces makes lattice simulations nearly impossible due to sign and phase oscillations. Given 
these difficulties we try a different approach. We keep the lattice spacing large, a ~ (100 
MeV) -1 ~ 2 fm, and again consider chiral effective field theory at leading order. But 
this time we introduce higher-derivative operators which improve the S-wave scattering 
amplitude at higher momenta. We expect that this should remove the clustering instability 
in the four-nucleon system. 
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FIG. 3: Suppression of finite cutoff errors by broadening the leading order contact interactions. 



X. HIGHER-DERIVATIVE TERMS 

As explained in the previous section, the interactions at leading order with delta function 
contact interactions are too strong at large momenta. Their contribution would be appro- 
priately weakened by interactions of higher chiral order. But a full investigation of higher 
order contributions is beyond the scope of this first exploratory study and is deferred to 
future work. Instead we consider here the effect of higher derivative terms which reduce 
cutoff errors by improving the delta function contact interaction. We fix the problem of 
clustering instability by introducing an (9(a 2 )-improved broadening for the leading contact 
interactions C and Cj, as illustrated in Fig. [31 This is by no means a full NLO calculation, 
but rather an LO calculation with 0(a 2 )-improvement to reduce cutoff errors. 
To this aim, we define the momentum- dependent densities, 




(64) 




(65) 
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where q s is the spatial momentum on the lattice. We can write q s as 



(66) 



where the components of k s are integers from to L — 1. 

The transfer matrix with only nucleon fields was defined in (H5l) . The contact 

interactions in M^ nt ' have the form 



7=1,2,3 n s 



q s L 7=1,2,3 

We replace these by the momentum-dependent interactions, 



(67) 



J/E^) -lca t p at ' a (q s )p a '' a (-q s ) - l -Cja t £ P f' a (q s )pf' a (-q s 

qs I 7=1,2,3 

where the coefficient function f(q s 2 ) is defined as 



(68) 



f(q s 2 ) = f 1 exp 



~ h E ( x - cos ^) 



/ s =l,2,3 

and the normalization factor fo is determined by the condition 



(69) 



/o-TlE 



cxp 



E ( x ~ cos ^) 

^=1,2,3 



(70) 



The coefficient b is determined at a later stage when we find the effective range. For small 
q s we see that f(q s 2 ) reduces to a Gaussian form, 



(71) 



We can introduce exactly the same momentum- dependent interactions in the transfer 
matrix formalism with auxiliary fields, 



Z LO oc J Dtt'jDsDs! exp [-3^ - S ss ] x Tr {Af (£t-1) (7^, s, sj) x • • • x M (0) (7Tj, s, sj)} . 
To do this we replace 



(72) 
(73) 



J n 
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by the nonlocal action 

- s^s^^f^ins-n^sin^nt) + Yl s i(^s,rh)f~ 1 (n s -n' s )s I (n' s ,n t ). (74) 



n s ,n'n t I n s ,n'n t 



The function / is defined as 



f-\n s - n' s ) = -L £ _|_ e -^ (^-^). (75) 



9s 

When the auxiliary fields are integrated out we recover the momentum-dependent interac- 
tions in (1681) . 



XI. THE TWO-NUCLEON SYSTEM REVISITED 



Using the new transfer matrix with momentum-dependent interactions we now revisit 
the two-nucleon system. Just as before we set C3 Sl and Cig to give the physical deuteron 
binding energy and physical 1 S scattering length. We also tune the coefficient b so that 
when Cs Sl and &s are determined, we also get the correct value for the average effective 
range § (Vj 50 + r s A . For a = (100 MeV)^ 1 and a t = (70 MeV)" 1 we find Ca 5l = -4.780 x 
10" 5 MeV~ 2 , Ci 5o = -3.414 x 10" 5 MeV" 2 , and b = 0.6. The new results are shown in 
Tables 3 and 4. 

Table 3: S-wave parameters 





rj 5 ° (fm) 


a s ?att (fm) 


r 3 Sl (fm) 


lattice 


3.20(1) 


5.30(1) 


1.46(3) 


experiment 


2.75(5) 


5.424(4) 


1.759(5) 



Table 4: Properties of the deuteron 





r d (fm) 


Q d (fm 2 ) 


lattice 


1.989(1) 


0.278(1) 


experiment 


1.9671(6) 


0.2859(3) 



The agreement with experimental values is now good. There is clear improvement over the 
earlier results shown in Tables 1 and 2. 

We can probe the shape of deuteron wavefunction by computing the nucleon density 
correlation function 

(:p at ' a (ryp^(0):). (76) 
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If A is the total number of nucleons then 

A = J2(p at > a (n s )). (77) 

We also find 

£ (■■ /'°K)/' a (0) :} = ^ (/' a K)/'°(0)} - (/' a (0)) 

fis fis 

= L- 3 (A 2 - A) . (78) 

Let us define the normalized density correlation function as 

G pp (n s ) = L 3 (A 2 - Ay 1 (: p at > a (n s )p a ' > a (0) :) . (79) 

In Figure HI we show G pp (n s ) for the deuteron in the xy-pl&ne. We have aligned the deuteron 
spin in the +z-direction. Keeping the deuteron spin in the +z-direction, in Figure Owe show 
Gp P (n s ) in the T/z-plane. A small asymmetry can be seen between the y and z directions. 
This is a signal of the deuteron quadrupole moment and can be seen more easily in Fig. [6l 
where we have taken an antisymmetric combination under interchange of y and z . 



XII. TRANSFER MATRIX PROJECTION METHOD FOR LIGHT NUCLEI 



We simulate light nuclei by using the Monte Carlo transfer matrix projection method 
introduced in 84( . Since this method may be unfamiliar, we first give an overview of the 
calculation using continuum notation before describing the details of the lattice transfer 
matrix calculation. 

Let |^/| e Ar) be a Slater determinant of free particle standing waves in a periodic cube for 
Z protons and N neutrons. We define A = Z + iV as the total number of nucleons. Let 
Hlq denote the Hamiltonian including instantaneous one-pion exchange and the improved 
higher- derivative contact interactions. Let HsutM De the same Hamiltonian but with both 
Gj and set to zero. As the notation suggests, Hsu(A)f is invariant under 577(4) Wigner 
symmetry. Wigner symmetry refers to an idealized limit where spin and isospin degrees of 
freedom can be interchanged and the 577(2) x SU (2) spin-isospin symmetry is elevated to 
an 577(4) symmetry. Let us define a trial wavefunction 

|¥(f)> = exp [-H su{ a)A (80) 
23 




2 3 4 _4-3 y (lattice units) 
x (lattice units) 



FIG. 4: Density correlations in the rcy-plane for a deuteron with spin in the +z-direction. 



G pp (0j,z) 




2 3 4 _4-3 z (lattice units) 
)> (lattice units) 



FIG. 5: Density correlations in the yz-plane for a deuteron with spin in the +z-direction. 
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G pp (O,y,z)-G pp (O,z,)0 




2 3 4 _4-3 z (lattice units) 
y (lattice units) 



FIG. 6: A linear combination of density correlations in the yz-plane that is antisymmetric under 
interchange of y and z. The deuteron spin points in the +z-direction. 

With this trial wavefunction we define the amplitude, 

Z(t) = {V(t')\e W [-H LO t}Mt')), (81) 
as well as the transient energy, 

E(t) = ~\lnZ(t)]. (82) 

In limit of large t we get 

lim E(t) = E , (83) 

t— *oo 

where E is the energy of the lowest eigenstate |^o) of Hlo with a nonzero inner product 
with | \& (£')). In order to compute the expectation value of some normal-ordered operator 
O we define 

Z (t) = (*(0I exp [-H LO t/2] O exp [-H LO t/2] \V(t')) . (84) 
The expectation value of O for |\P ) can be computed in the large t limit, 

Km^ = (* |0|*o>. (85) 

t^oo Z (t) 
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In this two-step approach we use exp [—Hsu(4)rft'] as an approximate inexpensive low-energy 
filter and exp [—H LO t] as an exact low-energy filter. The projection exp [— Hsu{i)^t'] is 
computationally inexpensive because the path integral for leading-order pionless effective 



field theory in the SU(4) limit is strictly positive for any even number of nucleons 85|, l86 |. 
Although there is no positivity theorem for odd numbers of nucleons, sign oscillations are 
also suppressed in odd systems because it is only one particle or one hole away from an even 
system with no sign oscillations. 

In the lattice transfer matrix formalism we construct \^(t')) using 

)> = Jds exp [S ss (s)} x Mf^\s) x • • • x M^/s) \¥™ N ) , (86) 

where t' = L to a t . M^^(s) is the same as M^faj, s, sj) except with = Cj = 0. We 
have omitted the dependence on tc'j and si since these completely decouple in the SU(4) 
symmetric theory. The amplitude Z(t) is constructed using 

Z(t) = J DitjDsDsi exp [S nw - S ss ] x | M (Lt * _1) «, s, Sj ) x • ■ ■ xM^^'j, s, S/ ) \^{t')} 

(87) 

where t = L t .a t . In this case the full leading-order transfer matrix is used. We compute 
Zo(t) by inserting the normal-ordered operator O in the middle of the string of transfer 
matrices, 

Z {t) = J Dir'jDsDs! exp [S n7T - S ss ) x (#(t')| M^" 1 ^, s, s z ) x • • • 

... x M^W^s^s^OM^I^^s^S!) x ■■■ x M (0) (7r;,s, S/ ) \^{t')) . 

(88) 

A schematic overview of the transfer matrix calculation is shown in Fig. [71 

Let j^i) , • • • , \ipA) be the free particle standing waves comprising |*|^r). We note that 
for the transfer matrices in the auxiliary field formalism there are no direct interactions 
between particles, only single-nucleon operators interacting with the background pion and 
auxiliary fields. It is easier to see this fact if we pretend for the moment that each of the 
A nucleons has an extra quantum number which makes them distinguishable. Let us label 
the extra quantum number as X, where X = 1, • • • , A. Then we have 



a i,j> a i,j ~^ \ a i,j,x, a i,j,x \ ■ (89) 
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operator insertion for 
expectation value 

<*| e w e | SU(4)/ fullLO (/ M1LO 5[/(4)/|^| e w e > 



2L t + L t . L, + L t . L t + LJ2 L t 

'o 'i l o n 'o n 'o 



FIG. 7: Overview of the various pieces of the transfer matrix calculation. 



We use an X subscript to indicate this replacement for creation and annihilation operators. 
The transfer matrices Mg^)^ and M*™*) factorize into transfer matrices for each X, 

II M svlmx> (9°) 

X=l,- ,A 
X=X,- ,A 

So long as the initial and final state wavefunctions are completely antisymmetric in X, 
this error in quantum statistics has no effect on the final amplitude. We can write the full 
A-nucleon transfer matrix element as a determinant of an A x A matrix of single-particle 
matrix elements, 

MM, s, sr, t>, t) = (A,x\ M^ L r 1] x • • • x Mf^M^ +L ^ 1] 

• ■ ■ x M^M^] X x ... x M® WtX \^,x) • (92) 

The indices i, j go from 1 to A. The calculation of Z(t) now reduces to computing 

Z(t) = J Dn'jDsDsj exp {-3^ - S ss } det M{ttj, s, s /; t', t). (93) 

There are many different ways to compute the amplitude Zo{t), and the most efficient 
method will depend on the operator O. Insertions of an operator measuring spatial difermion 
pair correlations was discussed in Here we consider an operator that measures spatial 
correlations of the total nucleon density, 

O = : p a '' a (n s )p at > a (Q) : . (94) 

For this operator it is convenient to use the identity 

0= lim — — M(e 1 ,e 2 ), (95) 
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where 

M(e 1 ,e 2 ) = : exp e 1 ^ a L'(^)Oij(n s ) + e 2 ^ a {j(Q) a i,j$) '■ • ( 96 ) 

hi i,j 

This is useful because M(ei,e2) itself looks like a transfer matrix with only single-nucleon 
operators. Let us define the new single-particle matrix elements, 

Mij(7r' I ,s,s I ,t',t,e 1 ,e 2 ) 

- {ipi,x\ M S u{A)j,x x • ■ ■ x M su(mx M x X • • • 

... x M^ /2) M x ( ei , e2 )M^ + ^ /2 - 1} x ... x M^M^] X x ... x mS (4V>x | fe > . 

(97) 

We then find 

d d F 

Z (t)= lim Dn'jDsDs! exp{-S n7T - S ss }detM(7i' I ,s,s I ,t l ,t,e 1 ,e 2 ). (98) 

€1,62— >Q <J€\ Ot 2 J 



We use hybrid Monte Carlo to update the fields it'j, s, sj [51]. We introduce conjugate 
fields P-k'^PsiVsj and use molecular dynamics trajectories to generate new configurations for 
Pn'^Ps, P SI , 7Tj, s, S! which keep 



h hmc = \Y,pl'S n ) + \ Erf(*o + \ E A w + v w> s ' ^ (") 

constant, where 

s, SJ ) = + S ss - log {|det M(n'j, s, sj, t', t end )\} . (100) 

tend denotes the largest value of t being considered. Upon completion of each molecular 
dynamics trajectory, we apply a Metropolis accept or reject step for the new configura- 
tion according to the probability distribution e~ BMC . This process of molecular dynam- 
ics trajectory and Metropolis step is repeated many times. Each time for the current 
configuration, C, we measure 

e ie (c) = detA^(7r^,g,g J ,t / ,t end ) 
|det.M(7r;, S , S/ ,t',t end )|' 

, _ det M(ir' T ,s,ST,t't) 

v ; |det.M(7r^s,sj,t',t erid )r v ; 



and 



ry , . detM(7r' I ,s,s I ,t',t end ,ei,e 2 ) f . 

Z ° (ei ' 62 ' C) " |det^K, S , S/ ,t',t end )| ■ (1 ° 3) 
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We take the ensemble averages for each measurement and form the ratios 

Z{t) (Z(t,C)) 
Z(t end ) (e*(C)) ' 

and 



(104) 



Zo(tena) _ Um d d (Z ( y e 1 ,e 2 ,C)) 



Z(t cnd ) e u e^o de 1 8e 2 (e^(C)> 
The ground state energy E is extracted using 



z{t , ^) =exp{EoAth (106) 

t ond ^oo Z (t end J 



and the expectation value (fy \ O \^ ) is calculated using 

,. Zo(t en d 
hm -— r 



l im Z o(tead) =( ^ Q | |^ o ^ (1Q7) 



XIII. NUMERICAL CHECKS USING THE TWO-NUCLEON SYSTEM 

We have already calculated properties of the deuteron using the transfer matrix with pions 
and auxiliary fields integrated out. In this section we make use of the two-nucleon system 
as a high-precision test of the transfer matrix Monte Carlo code. We calculate the same 
observables using both the Monte Carlo code and the exact transfer matrix, which is referred 
to as "Exact" method in the following. We use the lattice parameters defined previously, 
a = (100 MeV)" 1 , a t = (70 MeV)" 1 , C 3Sl = -4.780 x 10" 5 MeV" 2 , Ci So = -3.414 x 10~ 5 
MeV -2 , and b = 0.6. We set L = 3, L to = 2, L ti = 2. We consider a small system so that 
the stochastic error is small enough to detect disagreement at the 0.1% — 1% level. 

For the first test we consider |^/^ e ^r) with Z = 0, N = 2. The standing waves ipi,2 
comprising I^^jv) are 



i(n s ) \4>i) = L 3/2 ^ i0 5 i5 i, 
(0\a h3 (n s )\^ 2 ) = L-^ 2 5 t , 1 5 J , 1 . (108) 

This corresponds with a J z = 0, J = dineutron system with zero total momentum. We 
compute E(t) as well as the density correlation 

G pp (n s ) = L 3 (A 2 - A)' 1 (: /' a (n s )/'°(0) :) . (109) 
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From G pp (n s ) we can determine the root-mean-square radius trms as well as the quadrupole 
moment Q. The results for G pp (n s ) are shown in Table 5 and the results for E(t), trms, 
and Q are shown in Table 6. 

Table 5: G pp {n s ) for the J z = 0, J = dineutron system 



n s 


TV /T , /' l i 

Monte Carlo 


Exact 


(U, U, (J J 


0.09o0(o) 


0.09575 


(1,0,0) 


0.04506(5) 


0.04508 


(0,1,0) 


0.04515(6) 


0.04508 


(0,0,1) 


0.04500(5) 


0.04508 


(0,1,1) 


0.03354(4) 


0.03363 


(1,0,1) 


0.03365(4) 


0.03363 


(1,1,0) 


0.03358(6) 


0.03363 


(1,1,1) 


0.02876(4) 


0.02878 



Table 6: E(t), trms, and Q for the J z = 0, J = dineutron system 





Monte Carlo 


Exact 


E(t) (MeV) 


-5.84(9) 


-5.917 


r RM s (fm) 


1.402(4) 


1.4015 


|Q| (fm 2 ) 


< 10" 5 






The final values are not of physical relevance since the volume and number of time steps are 
very small. The important result is we find no disagreement between Monte Carlo results 
and the exact results beyond the stochastic error level. 

For the second test we consider I^IPat) with Z — 1, N — 1. This time the standing waves 
comprising |^^-) are 

(0|a M -(n s )|Vi) = L- 3/2 ^Ai, 

(0|a„(n,)|V 2 ) = ^ 3/2 ^Ao- (HO) 

This corresponds with a deuteron system with J z — 1, J = 1 and zero total momentum. 
The results for G pp {n s ) are shown in Table 7 and the results for E(t), trms, and Q are shown 
in Table 8. 
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Table 7: G pp (n s ) for the J z = 1, J = 1 deuteron system 



/ \ 

\Jl x , Thy, 7l z ) 


Monte Carlo 


Exact 


(0,0, 0) 


0.1230(3) 


0.12262 


(1,0,0) 


0.04233(4) 


0.04240 


(0,1,0) 


0.04247(5) 


0.04240 


(0,0,1) 


0.05563(5) 


0.05572 


(0,1,1) 


0.03415(5) 


0.03422 


(1,0,1) 


0.03423(3) 


0.03422 


(1,1,0) 


0.02764(4) 


0.02766 


(1,1,1) 


0.02650(3) 


0.02649 



Table 8: E(t), r RMS , and Q for the J z — 1, J — 1 deuteron system 





Monte Carlo 


Exact 


E(t) (MeV) 


-9.26(9) 


-9.311 


r RM s (fm) 


0.6957(2) 


0.69564 


Q (fm 2 ) 


0.1026(2) 


0.10283 



Again the final values are not of physical relevance since the volume and number of time 
steps are small. We find no disagreement between Monte Carlo results and the exact results 
beyond the stochastic error level. 

XIV. RESULTS FOR THE TRITON 

In our calculations we assume isospin symmetry and so our results for helium-3 and 
the triton will be the same. We choose to compare with experimental data for the triton 
since it is not affected by electrostatic repulsion. For our simulations we use the lattice 
parameters a = (100 MeV)" 1 , a t = (70 MeV)" 1 , Ca Sl = -4.780 x 10~ 5 MeV" 2 , G So = 
-3.414 x 10" 5 MeV~ 2 , and b = 0.6. We set L = 5, L to = 8 and vary L u from 2 to 10. The 
standing waves comprising j^l^r) are 

(0\aij(n s ) = L-V 2 6 ifi 6 jtl , 
(0\a hJ (n s )\^ 2 ) = L~ 3 / 2 5 h0 5 jfi , 

(0\a iJ (ii 8 )\fo) = L- 3 '\ 1 6 j , 1 . (Ill) 
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FIG. 8: The transient energy for the triton system as a function of imaginary time t. 

This corresponds to a triton system with J z — |, J = | and total momentum zero. 

For each value of L tl a total of about 5 x 10 6 hybrid Monte Carlo trajectories are gener- 
ated by 2048 processors, each running completely independent trajectories. Averages and 
stochastic errors are computed by comparing the results of all 2048 processors. In Fig. [8] 
we show the triton energy as a function of imaginary time t. The error bars denote the 
stochastic error. To remove the transient signal of higher energy states, we fit E(t) to a 
decaying exponential form at large t, 

E(t) -> E + ce- AEt . (112) 

A least squares fit gives an asymptotic value Eq = —8.9(2) MeV. Our result is within 
5% agreement with the experimental value of —8.48 MeV. Note that the triton appears to 
be overbound. However, for the triton the lattice volume of (9.85 fm) 3 that we use might 
still be too small and finite volume effects can possibly occur. In this case increasing the 
lattice volume would lead to a slightly less bound triton. Both finite volume effects and the 
systematic inclusion of higher chiral orders in the effective potentials will be investigated in 
future work. 

We also measure the nucleon density correlation G pp {n s ). In Fig. [9] we show G pp {n s ) 
in the x?/-plane at imaginary time t = 0.143 MeV -1 . From G pp (n s ) we can also extract 
the root-mean-square radius for the total nucleon density. In Fig. [10] we show the triton 
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FIG. 9: The nucleoli density correlation G pp (n s ) for the triton in the xy-plane at imaginary time 
t = 0.143 MeV- 1 . 

root-mean-square radius for the nucleon density as a function of imaginary time t. We 
again fit to a decaying exponential form, 

rRMs(t) -> r RMS + ce- Am , (113) 
and extract an asymptotic value of 2.27(7) fm. This is about 30% larger than the ex- 



perimental value of 1.755(9) fm for the root-mean-square radius of the electric charge 88] 



and the point proton root-mean-square radius 1.60 fm 46]. Similar values for the triton 



root- mean-square radius are also obtained in the pionless framework 89]. 



XV. RESULTS FOR HELIUM-4 



For our simulations for helium-4 we use the same lattice parameters a = (100 MeV) 
at = (70 MeV)- 1 , Cz Sl = -4.780 x 10~ 5 MeV" 2 , Ci So = -3.414 x 10~ 5 MeV" 2 , and b = 0.6. 
We again set L = 5, L to = 8 and vary L t . from 2 to 10. This time the standing waves 
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FIG. 10: The triton root-mean-square radius for total nucleon density versus imaginary time t. 



comprising l^f* 3 ^) are 

(0| aij(n s ) \ipx) 
(0| a i:j (n s ) \ij} 2 ) 
(0| a itj (n s ) |V> 3 ) 
(0| 0ij{n 8 ) \ip4) 



= L~ 3 / 2 5 ij0 8j t i, 
= L 3 ^ 2 5j i o5j,o, 
= L~ 3 ' 2 5^x5^1, 

= L-^ 2 8 iA 5j, . (114) 



This corresponds with a helium-4 system with J z = 0, J = and total momentum zero. 
For each value of L t . we again produce about 5 x 10 6 hybrid Monte Carlo trajectories using 
2048 processors running independent trajectories. 

In Fig. [TT] we show the energy for helium-4 as a function of imaginary time. If we fit 
E(t) to a decaying exponential form at large t we find an asymptotic value of —21.5(9) 
MeV. This is about 25% smaller in magnitude than the experimental result of —28.296 
MeV. Note that we compare directly with the experimental value and have not corrected 
for small corrections due to electromagnetic effects. In Fig. [T5] we show the nucleon density 
correlation G pp (n s ) in the xy-plane at imaginary time t = 0.143 MeV -1 . From G pp {n s ) 
we compute the root-mean-square radius for the total nucleon density, and in Fig. [13] we 
show the helium-4 root-mean-square radius for the total nucleon density as a function of 
imaginary time t. Fitting to a decaying exponential form we find an asymptotic value of 
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FIG. 11: The transient energy for helium-4 as a function of imaginary time t. 
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FIG. 12: The nucleon density correlation G pp (n s ) for helium-4 in the xy-plane at imaginary time 
t = 0.143 MeV" 1 . 
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FIG. 13: The helium-4 root-mean-square radius for total nucleon density versus imaginary time t. 



'"rms = 1.50(14) fm. This is within 10% of the experimentally observed value of 1.673(1) fm 
for the root-mean-square radius for the electric charge {9o| and the point proton root-mean- 



square radius 1.47 fm 



46|. 



XVI. DISCUSSION 



A. Room for more improvement 

The lattice calculations in this study used leading-order chiral effective field theory with 
improved contact interactions of the form 



l -Ca tP a "' a (q s )p a '' a (-q s ) - \Cja t £ pf > a (q s )pf > a (-q s 



7=1,2,3 



where 



Z s =l,2,3 



exp 



Is 



/ s =l,2,3 



cos q ls 



;ii5) 



(116) 



117) 



The three unknown constants C,Cf,b were determined using the deuteron binding energy, 
1 Sq scattering length, and the average effective range | (r So + r Sl ^j. The scattering lengths 
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and effective ranges were computed using Liischer's formula. We used only one value for 
the lattice spacing, a = (100 MeV)" 1 and a t = (70 MeV)^ 1 . The simulations for triton 
and helium-4 were done at only one volume, a cubical lattice with length L = 5 or 9.85 fm. 
In future studies both the dependence on the lattice spacings and finite volume should be 
investigated in some detail. 

Lattice results for the deuteron root-mean-square radius and quadrupole moment agree 
with experimental values at the 5% level. The lattice calculation of the triton binding energy 
agrees with experiment to within 5%, while the triton root-mean-square radius is larger by 
about 30%. For helium-4 we find the binding energy is about 25% smaller than experiment, 
while the root-mean-square radius agrees within 10%. This is clear improvement over the 
first attempt with zero-range contact interactions which led to a clustering instability in the 
four-nucleon system. However, in order to make a more rigorous statement we also have to 
consider the full set of interactions which arise at next-to-leading order in chiral effective 
field theory. We will discuss this procedure in future work. 

B. Computational scaling 

The lattice transfer matrix algorithm has several subroutines which scale differently with 
the number of nucleons A, the spatial volume L 3 , and the total number of time steps 
L t = 2L to + L ti . Multiplication of the transfer matrices at each time step on the one-particle 
states scales as Ax L 3 x L t . Constructing the Ax A matrix M(tc'j, s, sj,t', t) from the inner 
products of the one-particle states scales as A 2 x L 3 , while calculating the determinant of 
M^'j, s, si, t', t) scales as A 3 using LU decomposition. The molecular dynamics trajectories 
for the hybrid Monte Carlo updates involves computing the derivatives of H^mc with respect 
to Pn'jjPs, p sn n ii s > s i- The calculation of these derivatives require various loops which scale 
as L 3 x L t , L 6 x L t , and A 2 x L 3 x L t . 

With the same lattice parameters used in our triton and helium-4 simulations, we investi- 
gate the computational scaling for the transfer matrix algorithm as a function of the number 
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of nucleons A. We use L = 5, L to = 8, L ti = 10 and the free particle standing waves 



<o| 


a i,3 


(n s ) 


IV>1> 


= L 3/2 5i, 5j,i, 


<o| 


id 




|V>2> 




lr\ 1 

(o| 






1 / \ 


r-3/2r r 


(0| 


a i,j 




1^4) 


= L~ 3/2 5 ii i5 jfi , 


(0| 




(n s ) 


1^5) 


= L- 3 / 2 2 1 / 2 cos( 


(0| 


a i,j 


(n s ) 


|V>6> 


= L- 3 / 2 2 1 / 2 cos( 


(0| 


a i,3 




|V>7> 


= L- 3 / 2 2 1 / 2 cos( 


(0| 


a i,3 


(n s ) 


|V>8> 


= L- 3 / 2 2 1 / 2 cos( 



2nn 2 
2ixn, 



118) 



The initial state |^2 e Ar) is composed of all free particle standing waves |?/>j) with i < A. For 
A = 2 this is the deuteron system, for A = 3 the triton system, for A = 4 helium-4, for 
A = 5 helium-5, for 1 = 6 lithium-6, for A = 7 lithium-7, and for A = 8 beryllium-8. We 
note that for A > 4 the momentum of l^^ 6 ^) is not exactly zero but rather a wavepacket 
with a small spread in momentum centered around zero momentum. 

In Fig. [14] we show the CPU times for A = 2, • • ■ ,8 relative to the CPU time for A = 2 
with the same number of hybrid Monte Carlo trajectories. We see that for A < 8 the 
CPU time is approximately linear in A. This suggests that the subroutines which scale as 
A x L 3 x L t dominate the CPU time for the simulation for A < 8. 

In Fig. ITSlwe show the average phase (e lB ^ as defined in Eq. fllOip versus the number of 
nucleons. We note the relative maxima in (V 61 ) at multiples of 4. This can be explained by 
the suppression of sign and phase oscillations due to approximate 577(4) symmetry. The 
SU (4) suppression of oscillations is strongest for systems with equal numbers of spin-up and 
spin-down protons and neutrons. The results for the CPU time and average phase indicate 
that simulations of light nuclei with A < 8 can in fact be performed without much difficulty 
using the Monte Carlo transfer matrix algorithm presented here. We plan to pursue these 
studies in future work. 
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FIG. 14: CPU time versus the number of nucleons, A, measured relative to the A = 2 deuteron 
system. 
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FIG. 15: The average phase (e %e ) versus the number of nucleons, A. 



XVII. SUMMARY 



We have described simulations of light nuclei on the lattice using a transfer matrix pro- 
jection Monte Carlo method at leading order in chiral effective field theory. We included 
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lattice pion fields and auxiliary fields to reproduce the physics of instantaneous one-pion ex- 
change and the leading-order S-wave contact interactions. To avoid a clustering instability 
we also included higher-derivative contact interactions which adjust the S-wave scattering 
amplitude at higher momenta. There are a total of three unknown constants, C, Cj, b, 
and these were determined using the deuteron binding energy, 1 Sq scattering length, and 



experimental data at the 5% level for all calculated properties of the deuteron. The lattice 
result for the triton binding energy agrees with experiment to within 5% and the triton root- 
mean-square radius is within 30%. For helium-4 the binding energy is within 25% while the 
root-mean-square radius agrees within 10%. We expect that the description will improve 
when higher-derivative operators are treated systematically by matching phase shifts on the 
lattice at higher momentum. We have determined that the simulations for light nuclei with 
up to eight nucleons can be done without much difficulty using the lattice methods described 
here. 
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